function [F, G] = sim_moments(est, param)

    % add param
    define_param_est
    
    % constants
    define_omega_dist %firm type and density
    define_const %qgrid, phi_v, phy_y, D_F and other constants
    
    % solve the equilibrium
    eqm = solve_eqm_outer(param, const);    
    define_wage_schedule %add wage schedule

        % check
        if eqm.nq(param.Q)>0            
            disp('error: hit boundary')                                  
            [eqm.nq(param.Q)]
        end 
        scatter(const.Qgrid, eqm.nq, 'Marker', '.')
        xlabel('q')  
        drawnow limitrate      
    
    % calculate model moments   
    model = cal_model_moment(param, const, eqm); 
    M1 = model.mean_indegree_d;
    M2 = model.mean_outdegree_d;
    M3 = sqrt(model.var_lnsales_d);
    M4 = model.network_sales_share_d;
    M5 = model.exporter_share_d;
    M6 = model.export_intensity_d;
    M7 = model.avg_unwgt_lnwageS_d4;
    M8 = model.avg_wgt_lnwageS_d4;
    
    % return
    F = [M1,M2,M3,M4,M5,M6,M7,M8];
    G = [eqm.iter_inner, eqm.iter_outer];

end